package fem2.material;

import static java.lang.Math.log;
import static java.lang.Math.pow;
import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import fem2.Element;
import fem2.NotImplementedError;
import fem2.element.StructuralElement;
import fem2.enu.State;

public class NeoHookeanMM extends ElasticMaterialModel {

	public NeoHookeanMM(double E, double nu, double thickness, double rho, State ss) {
		super(E, nu, thickness, rho, ss);
	}

	@Override
	public IMatrix getConstitutiveMatrix(Element e, int k) {
		int dim = e.getMeshPart().getDimension();
		int ntrs = dim * (dim + 1) / 2;
		double[] xi = e.getIntegrationPoint(k);
		IVector epsilon = ((StructuralElement) e).computeStrain(xi);
		IMatrix C = Array1DMatrix.createRowMajor(ntrs, ntrs);
		computeC(C, epsilon);
		return C;
	}

	@Override
	public void Initialize(Element e) {
	}

	@Override
	public void InitializeSolutionStep(Element e) {
		/**
		 * Calculate first stress for the element
		 */
		StructuralElement se = (StructuralElement) e;
		for (int i = 0; i < se.getNumIntegrationPoints(); i++) {
			IVector sigma = se.getStress(i);
			double[] xi = se.getIntegrationPoint(i);
			IVector epsilon = se.computeStrain(xi);
			computeStress(sigma, epsilon);
		}
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		/*
		 * Calculate material response
		 */
		StructuralElement se = (StructuralElement) e;
		for (int i = 0; i < se.getNumIntegrationPoints(); i++) {
			IVector sigma = se.getStress(i);
			double[] xi = se.getIntegrationPoint(i);
			IVector epsilon = se.computeStrain(xi);
			computeStress(sigma, epsilon);
		}
	}

	@Override
	public void FinalizeSolutionStep(Element e) {
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

	private void computeStress(IVector s, IVector e) {
		switch (state) {
		case PLANE_STRAIN:
			double e11 = e.get(0);
			double e22 = e.get(1);
			double e12 = e.get(2);

			double s11 = ((lambda + 2 * e22 * lambda)
					* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) - 2 * mu
					* (-2 * e11 * (1 + 2 * e22) + pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1);
			double s22 = ((lambda + 2 * e11 * lambda)
					* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) - 2 * mu
					* (-2 * (e22 + 2 * e11 * e22) + pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1);
			double s12 = e12
					* (2 * mu - lambda * log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1);
			s.set(0, s11);
			s.set(1, s22);
			s.set(2, s12);
			return;
		case PLANE_STRESS:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		case THREE_DIMENSIONAL:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		default:
			throw new Error("invalid state");
		}
	}

	public void computeC(IMatrix c, IVector e) {
		switch (state) {
		case PLANE_STRAIN:
			double e11 = e.get(0);
			double e22 = e.get(1);
			double e12 = e.get(2);

			double c11 = -((4 + 8 * e22)
					* ((lambda + 2 * e22 * lambda)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) - 2 * mu
							* (-2 * e11 * (1 + 2 * e22) + pow(e12, 2))) * pow(2 + 4 * e22 + e11
					* (4 + 8 * e22) - 2 * pow(e12, 2), -2))
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (4 * (1 + 2 * e22) * mu + (2 + 4 * e22) * (lambda + 2 * e22 * lambda)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c12 = -((4 + 8 * e11)
					* ((lambda + 2 * e22 * lambda)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) - 2 * mu
							* (-2 * e11 * (1 + 2 * e22) + pow(e12, 2))) * pow(2 + 4 * e22 + e11
					* (4 + 8 * e22) - 2 * pow(e12, 2), -2))
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (8 * e11 * mu + 2 * lambda
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + (2 + 4 * e11)
							* (lambda + 2 * e22 * lambda)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c13 = 4
					* e12
					* ((lambda + 2 * e22 * lambda)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) - 2 * mu
							* (-2 * e11 * (1 + 2 * e22) + pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2)
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (-4 * e12 * mu - 2 * e12 * (lambda + 2 * e22 * lambda)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c22 = -((4 + 8 * e11)
					* ((lambda + 2 * e11 * lambda)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) - 2 * mu
							* (-2 * (e22 + 2 * e11 * e22) + pow(e12, 2))) * pow(2 + 4 * e22 + e11
					* (4 + 8 * e22) - 2 * pow(e12, 2), -2))
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (4 * (1 + 2 * e11) * mu + (2 + 4 * e11) * (lambda + 2 * e11 * lambda)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c23 = 4
					* e12
					* ((lambda + 2 * e11 * lambda)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) - 2 * mu
							* (-2 * (e22 + 2 * e11 * e22) + pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2)
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (-4 * e12 * mu - 2 * e12 * (lambda + 2 * e11 * lambda)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c33 = 4
					* (2 * mu - lambda * log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)))
					* pow(e12, 2) * pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2)
					+ (2 * mu - lambda * log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1) + 2 * lambda
					* pow(e12, 2) * pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1);

			c.set(0, 0, c11);
			c.set(0, 1, c12);
			c.set(0, 2, c13);
			c.set(1, 0, c12);
			c.set(1, 1, c22);
			c.set(1, 2, c23);
			c.set(2, 0, c13);
			c.set(2, 1, c23);
			c.set(2, 2, c33);
			break;
		case PLANE_STRESS:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		case THREE_DIMENSIONAL:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		default:
			throw new Error("invalid state");
		}
	}

	@Override
	public double computeEnergyAt(IVector e) {
		switch (state) {
		case PLANE_STRAIN:
			double e11 = e.get(0);
			double e22 = e.get(1);
			double e12 = e.get(2);

			/*
			 * TODO: check the equation below
			 */
			return (e11 + e22) * mu
					- (mu * log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2))) / 2.
					+ (lambda * pow(log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)), 2))
					/ 8.;
		case PLANE_STRESS:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		case THREE_DIMENSIONAL:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		default:
			throw new Error("invalid state");
		}
	}

	@Override
	public double computeEnergyAt(IVector e, IVector s) {
		/*
		 * TODO: check if there is better equation for this
		 */
		return computeEnergyAt(e);
	}
}
